BetheSalpeter equation in ABINIT¶
¶
These notes provide a brief introduction to the BetheSalpeter (BS) formalism outlining the most important equations involved in the theory. The approach used to compute the BS kernel and the macroscopic dielectric in an implementation based on planewaves and normconserving pseudopotentials is also discussed.
The conventions used in the equations are explained in the Manybody theory notes. A much more consistent discussion of the theoretical aspects of the BetheSalpeter equation can be found in [Onida2002] and references therein.
1 Optical properties and local field effects¶
Before discussing the BetheSalpeter problem, it is worth reminding some basic results concerning the ab initio description of optical properties.
In frequency and reciprocal space, the microscopic dielectric function is related to the irreducible polarizability by the following relation
from which the inverse dielectric function is obtained via matrix inversion.
The macroscopic dielectric function, \ee_M^{\text{LF}}(\ww), can be directly related to the inverse of a single element, the first (\GG_1 = 0, \GG_2 =0) element of the inverse of the microscopic dielectric matrix by means of:
The microscopic dielectric matrix is the one usually calculated within the RPA in a GW calculation. The optical absorption spectrum is simply given by the imaginary part of \ee_M^{\text{LF}}(\ww). Sometimes, especially when comparing with experimental spectra, the real part is simply called \ee_1 and the imaginary part \ee_2.
Note that above equation differs from
due to the so called localfield effects introduced by the presence of the crystalline environment. The former quantity with local fields (LF) is more accurate than the latter one without local field (NLF).
Important
The reason the two equations are different is because in the first, the expression in the denominator is the first element of the inverse of the whole microscopic dielectric matrix. This element depends on the entire matrix and cannot simply be calculated by taking the inverse of the first element of the noninverted matrix.
In the GW_notes, we have discussed how to calculate the irreducible polarizability and thus the absorption spectrum within the random phase approximation (RPA). It turns out, however, that the RPA dielectric function evaluated with KohnSham orbitals and eigenvalues yields absorption spectra that are in quantitative, and sometimes even in qualitative, disagreement with experiments. This is not surprising since the polarizability of the KohnSham system is not expected to reproduce the response of the fully interacting system.
Important discrepancies with experiment are found even when the DFT gap is corrected by replacing the KS eigenvalues by quasiparticle energies calculated in the GW approximation. This indicates that, as far as optical properties are concerned, there is some important physics that is not correctly captured by the GW approximation.
The fundamental reason for the failure of the RPA can be traced back to the neglect of the electronhole interaction that, in the manybody language, should be included through the vertex function. Replacing the vertex function with a local and instantaneous function is a too crude and unrealistic approximation for the manybody polarizability.
In the next section we will discuss how to obtain an improved approximation for the vertex and therefore an improved approximation for the polarizability that takes into account many body effects due to the electronhole interaction
2 The BetheSalpeter equation in a nutshell¶
A schematic picture of the different levels of the description of optical spectra available in ABINIT is given below.
The BetheSalpeter theory is formulated in terms of twoparticle propagators. These are fourpoint functions describing the motion through the system of two particles at once. We will label them L^0 and L, where the subscript zero indicates that the particles are noninteracting. By restricting the starting and ending point of the particles, the twopoint contraction of L^0 gives the reducible independentparticle polarizability according to
$$ \chi^0(12) = L^0(11, 22) $$ while the twopoint contraction of L gives the reducible manybody polarizability $$ \hat\chi(12) = L(11, 22) $$ that should not be confused with the irreducible polarizability of the manybody system (denoted with a tilde in the formulas in the GW_notes).
The utility of working in terms of the reducible quantities is that the macroscopic dielectric function with local field effects is obtained directly from the reducible polarizability using
Computing the reducible L directly, if possible, thus allows one to avoid the costly matrix inversion of the dielectric function that should otherwise be performed for each frequency.
One can show that L satisfies the Dysonlike equation:
that involves the BetheSalpeter kernel K:
where the bar symbol signifies that the Coulomb interaction has to be taken without its longrange Fourier component at \GG =0:
As discussed in [Onida2002], local field effects (LF) are included in the formalism through the exchange term while the Coulomb term describes the screened electronhole interaction that is responsible for the creation of the excitons. All interaction is in principle already contained in W, and the rewriting of the kernel in this way is mostly an algebraic trick which allows one to easily separate the contributions and calculate the optical spectrum both with and without the LF.
The Dysonlike equation for L becomes a matrix problem once a particular basis set is chosen to expand the fourpoint functions involved in the problem. In the standard approach, the polarisabilities and the BS kernel are expressed in the socalled transition space (i.e. as products of singleparticle orbitals) using:
and
where n_i is a shorthand notation to denote band, kpoint and spin index. The expansion is exact provided that the set of singleparticle orbitals from a complete basis set in the Hilbert space.
The main advantage of using the transition space for solving the problem is that the RPA polarizability is diagonal in this representation
thus leading to a considerable simplification in the mathematical formalism.
After some algebra (see [Onida2002] for a more complete derivation) one finds that, in a system with an energy gap, the Dysonlike equation for L can be expressed in terms of an effective twoparticle Hamiltonian, H , according to
The explicit form for H and F in the transition space is given by
where valence states are indicated by the indices v , v' while c , c' are used for conduction bands. The R submatrix is Hermitian and is usually referred to as the resonant block while C (the socalled coupling block) is a complex symmetric submatrix. Note that, due to the presence of the coupling term, the BS Hamiltonian is nonHermitian although the operator possesses a real spectrum.
The inclusion of spin in the formalism would require an extensive discussion on its own. Here, for simplicity, we limit the analysis to the case of spin unpolarized semiconductors (nsppol=1). In this particular case, the matrix elements of the resonant block are given by
with single particle transition energies on the diagonal.
The matrix elements of v and W are defined as:
Note, in particular, that only the static limit (\omega = 0) of W is involved in the last expression.
The coupling matrix elements are usually smaller than the resonant ones. This is especially true in bulk systems due to the different spatial behavior of conduction and valence states. In solid state calculations, it is therefore common practice to ignore the C block (the socalled TammDancoff approximation). Under this assumption the excitonic Hamiltonian is given by a Hermitian operator.
The variable bs_coupling defines whether the coupling block should be included in the calculation.
The macroscopic dielectric function is obtained by contracting the manybody L and then taking the optical limit of the \GG = \GG' = 0 component along a particular direction in \qqspace. The final result reads:
where
is the matrix element of the dipole operator in transition space.
By default the code calculates the macroscopic dielectric function taking the limit along six different directions in q space (the three basis vectors of the reciprocal lattice and the three Cartesian axis). It is possible to change the default directions using the variables gw_nqlwl and gw_qlwl.
3 Solving the BetheSalpeter problem¶
At this point it is clear that the evaluation of the macroscopic dielectric function within the BS formalism is a twostep process:

Construction of the H matrix in the transition space.

Computation of the macroscopic dielectric function using the two equations reported at the end of the previous section.
The flowchart of a typical BetheSalpeter run is schematically depicted in the diagram below:
The WFK file is represented with an ellipsis. The path on the left indicated with blue arrows represents the RPA calculation (optdriver=3) that produces the SCR file (see also the GW1 tutorial). Once the WFK and the SCR file are available, we can finally construct the Hamiltonian and solve the BetheSalpeter problem (the rectangle at the bottom of the flowchart).
For BS computations, it is common practice to simulate the selfenergy corrections by employing the scissors operator whose value can be obtained either from experiments or from abinitio calculations. The scissors operator allows one to avoid a costly GW calculation that should performed for all the \kkpoints and bands included in the transition space (the optional path on the right indicated with yellow arrows that corresponds to optdriver=4).
The construction of the Hamiltonian matrix represents a significant portion of the overall CPU time due to the large number of transitions needed for an accurate description of the frequencydependence of the polarizability. On the other hand, also the calculation of the macroscopic dielectric function poses numerical difficulties since an expensive matrix inversion should be performed for each frequency.
The code implements different methods proposed in the literature to avoid the matrix inversion for each frequency. The variable bs_algorithm is used to select among the different possibilities.
bs_algorithm=1 employs standard (Sca)LAPACK routines to obtain the spectral representation of H in terms of eigenvalues and right eigenvectors of H:
Then, as discussed in [Albrecht1998], the inverse of H  \ww is obtained according to
We do not elaborate this point further since the direct diagonalization is advantageous only in particular circumstances, as explained in the documentation.
bs_algorithm = 2 avoids the diagonalization using an iterative algorithm that constructs a basis set in which H is represented by a real symmetric tridiagonal matrix [Haydock1980]. Without entering into detail, one can schematically represent the Haydock technique as an algorithmic procedure that transforms a dense (hermitian) matrix into a sparse (tridiagonal) one:
Once the coefficient of the tridiagonal form are know, the macroscopic dielectric function is evaluated in terms of the continued fraction:
where a small complex shift, defined by zcut, is used to avoid the singularities along the real axis. The number of iteration required to converge is almost independent on the size of the problem, usually of the order of 100200 iterations. The algorithm involves simple matrixvector multiplications that are efficiently performed with BLAS calls and, most importantly, are easy to parallelize. Another distinct advantage is that the Haydock method is less memory demanding than the direct diagonalization since only three temporary vectors are required instead of the full set of eigenvectors.
Note that the original formulation of the method presented in [Haydock1980] assumes an Hermitian operator thus it can be used only for TDA calculations. We refer the reader to [Gruning2009] for a generalization of the method to the case in which the coupling term cannot be neglected. The main drawback of the method is that it does not give direct access to the excitonic spectrum hence it cannot be used to calculate binding energies or to plot the excitonic wavefunctions. Moreover, for the time being, bs_algorithm = 2 cannot be used for calculations in which the coupling term is included.
bs_algorithm = 3 employs a standard conjugategradient method to calculate the lowest eigenstates of the Hamiltonian. At present, this method is still under testing and does not support calculations with the coupling term.
4 Kernel matrix elements in reciprocal space¶
Our implementation employs planewaves to expand the periodic part of the Bloch states, u , and the twopoint function W(\rr,\rr') that becomes a \qqdependent matrix in reciprocalspace. The conventions used for the transforms are documented in this section of the GW notes.
The matrix elements of the exchange term are evaluated in reciprocal space using:
while the Coulomb term is calculated according to
The number of \GGvectors in W and in the modified Coulomb interaction is specified through ecuteps while the wavefuctions entering the oscillator matrix elements are expanded on a \GGsphere of energy ecutwfn. The computation of the oscillator matrix elements is discussed in this section of the GW Notes.
The input variable bs_exchange_term can be used to disable the computation of the exchange term, this option is mainly used for performing optical calculations without local field effects. The variable bs_coulomb_term governs the computation of the Coulomb term, the most CPUdemanding part due to the presence of the double sum over \GG vectors.
It is also important to stress that, in the two equations above, the kpoint index runs over the full Brillouin zone hence the size of the Hamiltonian is defined by the number of point in the full Brillouin zone and not by the number of points in the irreducible wedge.
5 Matrix elements of the dipole operator¶
The accurate calculation of optical properties require the correct treatment of the optical limit (\GG = 0, \qq \rightarrow 0) of the oscillator matrix elements. The computation of these terms deserves some clarification, due to the presence of the fully nonlocal pseudopotential term in the KohnSham Hamiltonian.
A linear expansion up to the first order in q of the oscillator matrix element results in:
where we have assumed b_1 \neq b_2 and the difference between the two wave functions at \kk  \qq and \kk has been neglected because it only introduces terms that are quadratic in \qq.
Unfortunately, the above expression cannot be directly used because the matrix elements of the position operator are illdefined when wavefunctions obey Bornvon Karman periodic boundary conditions. For this reason, the first order contribution has to be evaluated using the equivalent expression [Baroni1986]
The term involving the momentum operator is efficiently evaluated in reciprocal space with linear scaling in %npwwfn while the treatment of the nonlocal part of the pseudopotential is more involved and much more CPUdemanding.
The role played by this additional term is usually marginal in the case of GW calculations: the QP corrections are obtained by performing an integration in \qqspace and only the \qq \rightarrow 0 component of the inverse dielectric matrix is affected by the commutator of the nonlocal part of the pseudopotential.
For this reason it is common practice, especially during the GW convergence tests, to neglect this contribution by setting inclvkb = 0. Strictly speaking, however, this approximation is justified only in the case of calculations in bulk crystals provided that the BZ sampling is well converged. Particular care has to be taken when performing GW calculations in nonperiodic systems due to the reduced dimensionality of the BZ.
Important
The commutator of the nonlocal part should ALWAYS be included when studying optical properties, both at the RPA and at the BS level.